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Abstract. We show the existence of gap-Townes solitons for the multi- 
dimensional Gross-Pitaeviskii equation with attractive interactions and 
in two- and three-dimensional optical lattices. In absence of the peri- 
odic potential the solution reduces to the known Townes solitons of the 
multi-dimensional nonlinear Schrodinger equation, sharing with these 
the propriety of being unstable against small norm (number of atoms) 
variations. We show that in the presence of the optical lattice the so- 
lution separates stable localized solutions (gap-solitons) from decaying 
ones, characterizing the delocalizing transition occurring in the multidi- 
mensional case. The link between these higher dimensional solutions and 
the ones of one dimensional nonlinear Schrodinger equation with higher 
order nonlinearities is also discussed. 



1 Introduction 

One interesting phenomenon occurring in ultracold atomic gases trapped in pe- 
riodic potentials is the possibility to localize matter in states which can stay for 
a long time due to an interplay between nonlinearity, dispersion and periodic- 
ity. Such states (also called gap-solitons) have been observed in Bose-Einstein 
condensates (BEC) and in arrays of nonlinear optical waveguides [1121314) . For 
attractive atomic interactions in BEC and in absence of a periodic potential, sta- 
ble localized solutions are possible only in a one-dimensional (ID) setting since 
in two (2D) and three (3D) dimensions the phenomenon of collapse appears [5]. 
More precisely, one observes that when the number of atoms exceeds a critical 
threshold, the solution collapses in a finite time (blow-up) while for number of 
atoms below the critical threshold there is an irreversible decay of the state into 
background radiation. For the higher dimensional nonlinear Schodinger (NLS) 
equation, however, it is known that there exists an unstable localized solution, 
the so called Townes soliton [5] , which separates decaying solutions from collaps- 
ing ones. Townes soliton, however, exists only for a single value of the number 
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of atoms, being unstable against fluctuations around it (for slightly overcritical 
or undercritical number of atoms the solution collapses or decays, respectively). 
The situation is drastically changed in presence of an optical lattice (OL). To 
this regard, it has been shown that stable 2D and 3D solitons can exists in OLs 
both in BEC and nonlinear optics contexts |7l8l9ll0fTTj . Moreover, it is known 
that while the periodic potential can only marginally shift the critical value for 
collapse, it can substantially move the delocalizing transition curve, thereby in- 
creasing the soliton existence range in parameter space from a single point to 
a whole interval |12j . The typical situation with 2D and 3D BEC solitons in 
OLs is therefore the following: in the parameter space the stable localized so- 
lutions are confined from above by the collapse curve and from below by the 
delocalizing transition curve, thus, in contrast with the one dimensional case 
where there are no limits for the existence of localized states, strict limitations 
for soliton existence appear in multidimensional cases. From this point of view 
it is clear that for possible experimental observation of multidimensional BEC 
solitons the parameter design becomes very important. Since the collapse curve 
is only marginally affected by the periodic potential, to enlarge existence ranges 
of solitons it is of interest to give a full characterization of the delocalizing curve 
in parameter space. 

The aim of this paper is just devoted to this, i.e. we characterize 2D and 
3D delocalizing curves of gap-solitons in terms of an unstable solution of the 
multidimensional Gross-Pitaeviskii equation (GPE), which we call gap-Townes 
soliton. This solution can be viewed as a separatrix (it separates gap soliton 
states from extended (Bloch) states) and reduces to the known Townes soliton 
when the strength of the OL goes to zero. Similar solutions were found also for 
the ID NLS equation with higher order nonlinearities in (13j , where they were 
called gap-Townes solitons, and in [14] where they were termed Townes solitons. 
Conditions for the occurrence of the delocalizing transition phenomenon of one- 
dimcnsional localized modes of several nonlinear continuous periodic and discrete 
systems of the nonlinear Schrodinger type were also recently discussed in [15j . For 
the periodic multidimensional GPE the delocalizing curve has been characterized 
in [T2] as the critical threshold for the existence of one bound state in an effective 
potential. The characterization given here, however, is more general since it 
is valid also for ID NLS with higher order nonlinearities. To this regard we 
remark that in absence of confining potential the 2D and 3D GPE behaves 
similarly to the ID NLS with quintic and septic nonlinearities, respectively. The 
interplay between dimensionality and nonlinearity has been used to investigate 
collapse in lower dimensional NLS on the basis of pure dimensional arguments. 
In particular, the critical condition for collapse has been characterized as D{n — 
1) — 4 = 0, where n is the order of the nonlinearity in the equation and D 
is the dimensionality of the system [T6|. In the following we take advantage of 
this interplay to construct approximate gap-Townes soliton solutions of the GPE 
with multidimensional separable OLs, in terms of products of exact gap-Townes 
solutions of the ID NLS with higher nonlinearities. Remarkably, we find that, 
except for strengths of the optical lattices very small, our approach produces very 
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accurate gap-Townes solutions of multidimensional GPE with OL, thus giving 
an evident computational advantages. The results obtained in this paper can be 
seen as a generalization of the existence of gap-Townes solitons in the quintic 
NLS discussed in |13ll4j to the case of the multidimensional Gross-Pitaeviskii 
equation. 

We finally remark that the obtained results can also be applicable for pho- 
tonic lattices with Kerr type of optical nonlincarity where the existence of a 
critical threshold for the lattice solitons has been observed [17] . 

The paper is organized as follows. In section II we introduce the model equa- 
tions and discuss the link between multidimensional GPE with a separable trap- 
ping potential and the corresponding ID NLS equation with higher order non- 
linearity. We use a self consistent approach to approximate gap-Townes solitons 
of the GPE with products of exact gap-Townes soliton of the corresponding ID 
NLS equation with higher order nonlinearity. In Section III we discuss the ex- 
istence of localized solutions in the multidimensional GPE with OL by means 
by a variational approach and compare 2D and 3D results with those obtained 
from the VA applied to the quintic and septic NLS, respectively. In section IV 
we perform a numerical investigation of the existence (delocalizing) threshold 
for gap-Townes solitons of the 2D and 3D GPE. Finally, in the last section we 
briefly summarize our main results. 

2 Model equations and existence of gap-Townes solitons 

Let us consider the following Gross-Pitaevskii equation in d-dimensions (d — 
1, 2, 3,) as a model for a BEC in an optical lattice 12J 



where denotes the d-dimensional laplacian, 2^ =1 cos(2irj) denotes a square 
optical lattice with strength e, 7 is the coefficient of nonlinearity, and Xi = x,y,z 
for z = 1,2, 3, respectively. Here we will be mainly interested in cases d = 2 and 
d = 3. The existence of localized solutions of the multidimensional GPE with 
periodic potential and positive and negative nonlinearities (atomic scattering 
lengths), has been previously investigated both by variational analysis and by 
direct numerical simulations. In the following we concentrate on a topic which 
was not discussed in previous works, namely the existence of gap-Townes solitons 
in the multidimensional GPE and its link to the phenomenon of delocalizing 
transition. Due to the instability properties of these solutions it is difficult to 
find them without an analytical guide. To this regard we take advantage of the 
fact that the periodic potential is separable and in spite of the nonlinearity of 
the system we look for factorized stationary solutions of the form 
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In 2D case (d = 2) the substitution of the factorized ansatz into Eq. ([T]) gives: 
^ + ^ + e[cos(2:z) + cos(2zj)] + 7 |0i| 2 |0 2 | 2 = (3) 

01 02 

This equation can also be written as 

+ecos(2x)0i + ^|0 2 | 2 |0i| 2 0i = -Mi0i, 

02 yiy + ecos(2y)0 2 + ^|0i| 2 |02| 2 2 =-M202, (4) 

with /ii = /i2 = £t/2. By assuming 0i = 02 = and adopting a diagonal 
coordinate x = y = £, we have that Eqs. (j4|) become equivalent to the following 
ID eigenvalue problem 

55 + ecos(2O0+||0| 4 = -f0- (5) 

From this we see that there is a link between the 2D cubic NLS and the ID 
quintic NLS which implies a rescaling of parameters as: 7 — > 5r, fj, — » ^. 

The above equations can be easily extended to the 3D GPE with periodic 
potential. In this case Eq. ([5]) will be replaced by the following ID NLS equation 
with septic nonlinearity 

a + £ cos(2O0+ J0| 6 = -|0, (6) 
from which we see that in this case parameters must be rescaled according to: 

' 3 ' 3 \ 

It is appropriate to mention that a factorized solution of the form with 
the components solutions of the nonlinear eigenvalue problem ([5]) cannot be 
an exact solution of the 2D or 3D GPE, since, due to the nonlinearity, the 
problem is obviously not exactly separable. On the other hand, by imposing the 
coincidence of the solutions along the diagonal axis may be a constraint for a 
reasonable approximate solutions of the 2D and 3D problems, especially when 
the nonlinearity is small. In analogy with the ID NLS with quintic nonlinearity 
investigated in Ref. |13j we expect that the delocalizing curve coincides with 
the existence curve of gap-Townes solitons for which it was shown that the 
critical number of atoms decreases with increasing the strength of the OL. This 
means that in a deep OL the effective nonlinearity required for the existence of 
a gap-Townes soliton is smaller and the problem may become effectively close to 
separable. 

To check the correctness of this argument we construct factorized solutions 
|2]) of the 2D GPE by means of a self-consistent method which allows to solve 
the ID quintic NLS eigenvalue problem exactly (similar results can be obtained 
for the 3D case). 

In the top panel of Fig. 1 we show a 2D gap-Townes soliton obtained from 
the product ansatz using exact (self-consistent) gap-Townes solitons of the cor- 
responding quintic NLS equation in Eq. (|5|). The lower left and right panels 
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Fig. 1. Gap-Townes soliton (top panel), gap soliton (lower left panel) and Bloch 
state (lower right panel) of the GPE for parameter values e — 5, 7 = 1. 
The critical number of atoms in normalized units for the gap-Townes soliton 
is N c = 0.4261. The gap soliton and the Bloch state in the lower panel are ob- 
tained for slightly overcritical (N=0.4347) and undercritical (N=0.4176) values, 
respectively. 
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show, respectively, the gap soliton and the extended Bloch state found at en- 
ergy slightly below and slightly above (bottom of the lowest band) the one of 
the gap-Townes soliton. To check the reliability of the factorized ansatz we have 
computed the time evolution under the original 2D GPE equation using the 
factorized solution as initial condition. This is shown in Fig. 2 where the time 
evolution of the gap-Townes soliton in Fig. 1 (central panel) and the ones ob- 
tained for slightly overcritical and undercritical numbers of atoms are shown. We 
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Fig. 2. Time evolution of the gap-Townes soliton shown in Fig. 1 with critical 
number of atoms (central panel) and for slightly overcritical (left panel) and 
undercritical (right panel) number of atoms. 



see that while the product solution constructed from the quintic NLS remains 
localized for a long time, a slight increase or decrease of the number of atoms 
produces shrinking or decay of the solution, respectively. This clearly shows the 
existence and role of gap-Townes solitons of the 2D GPE with optical lattice in 
characterizing the delocalizing threshold. 



3 Variational analysis and existence of localized states 



The existence of localized states in multidimensional GPE with OL and in ID 
NLS with higher order nonlinearities can be investigated by means of the vari- 
ational approach (VA). To this regard we first consider Eq. ^ and search for 
stationary solutions of the form i/j = Ue~ %tit and consider the Lagrangian density 
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associated with the equation for the stationary field U 

d 



^ cos(2a;i) 



i=l 



u 2 - lu\ 



By taking a Gaussian ansatz for U 



U = Ae 



(7) 



(8) 



and performing spatial integration we obtain the following effective lagrangian 
L e ff = J Cdx for parameters A, a 



L 



eff 



A*_ 

2 



(- 

V a 



H - dee~ 1/a 



A 2 j 
2i +1 



1,2,3. 



(9) 



Variational parameters A, a, indicating the amplitude and inverse width of the 
localized state, are linked to the number of atoms by the relation N = A 2 {^) d / 2 . 
From the conditions of stationarity of the effective lagrangian dL e ff/da = 
dL e ff/dA = 0, we get the following equations relating N with a and the chem- 
ical potential fi 



u= -a — dee x l a 
p 2 



7 iV (a\ d / 2 

~2~d/2 



(10) 



N 



2e 



-l/a 



d= 1,2,3. 



In Fig. [3] we depict the (N, //) curves obtained from the above transcendental 
equations for fixed values of e and for the cases d = 1,2, 3. We see that for the ID 
case dN/dfj, is always negative, this means, according to the Vakhitov-Kolokolov 
(V-K) criterion [18], that the solution exists and is stable for any value of e 
without limitations on the effective nonlinearity Nj. In the 2D threshold 
in N-f appears which is predicted by VA to be exactly Air for e — 0. Also notice 
that in this case dN/dfi < is still satisfied for most branch curves, meaning 
that the solution is usually stable. The situation is quite different in the 3D case, 
where there is no limiting threshold for existence but most of the curves display 
a positive slope meaning that the solution is unstable. In particular, from the 
upper inset of the right panel we see that for curves close to the \i — axis, 
dN/dfj, change signs at iV7 w 70 and becomes positive for higher values of Nj. 
The almost horizontal curves for lower values of -/V7 are displayed in the lower 
inset of the figure, from which we see that dN/dfi changes from negative to 
positive after the curves have reached a maximum at values of /1 which depend 
on e. 

Notice that for the 3D case the above equations predict for e = the number 
of atoms dependence on the chemical potential as 
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Fig. 3. Left panel. VA existence curves for the ID case (d = 1 in Eqs. US])) 
for different values of e. The inset shows an enlargement for small values of \x. 
Curves refer to e values increased in steps of 0.5 starting from e = 0.5 (upper 
curve and most left curve in the inset) to e = 2.0 (lower curve and most right 
curve in the inset). Central panel. Same as for the left panel but for the 2D case 
(d = 2 in Eqs. (fT5|) ). The inset shows an enlargement for small values of /x close 
to the threshold N = 4tt. Right panel. Same as for the left panel but for the 3D 
case (d = 3 in Eqs. (I15p ). The upper inset shows an enlargement of the curves 
close to the /x = axis, while the lower inset displays the curves close to the 
A^7 = axis, on a larger scale. 



Similar dependence N vs. fi was previously obtained for stationary solutions of 
3D NLS with cubic and quintic nonlincarity [19 . Also notice that according to 
V-K criterion this solution is unstable. From the condition N > 0, we obtain the 
limitation on the soliton width as 



- 2 z- 1/a < 1. 



(11) 



while for s >> 1 we obtain that the soliton exists if the width satisfies the 
condition 

2(log(<5) + log(log(£))) 

and a > a c2 ~ \/2e— 1/2. For e = 5, for example, we obtain a c \ — 0.2, while the 
exact value is a e » 0.17. 

In view of the analogy of the 2D and 3D GPE with the ID NLS with quintic 
and septic nonlinearity, respectively, it is of interest to compare the above VA 
equations with these cases. To this regard, we consider the ID GPE with a high 
order nonlinearity of the type 7|V ; | Q V' with a = 2,4,6 



iipt + ipxx + l\i>\ a i^ + £ cos(2a;)V' = 0. 



(13) 



Using the same approach as before, one can readily show that the effective 
lagrangian in this case is 



L ~ A ( 

L eff - — I ~ 



1/2 
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ee 
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3/2 



(14) 
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from which the following VA equations are derived 

° Ee -i/q ^ Na/2 rr /4 as) 

££ (£±2)1/2 UJ ' (15j 



a- 7 



here a = 2, 4, 6. We see that the case a = 2 coincides with the case d = 1 
considered above, and the case a — 4 with the quintic VA equations derived in 

Notice that for e = Eq. (fTB")) admits exact solutions also for a — 4, 6. For 
the septic case, indeed, we have, using ip = ue mit , m — — // > 0, that 

i> = ( — ) 1/6 sech 1/3 (3^x), (16) 
7 

is an exact solution with a norm 

_/V = — - {3> (17) 

3r(f) 7 V3 m i/6' ^'J 

where ^(x) is the gamma function. The Hamiltonian for this solution is equal 
to 

2 7, ,,8^_n,X /6 



ff= / (IV'.I 2 - J^| 8 )^ = 0.44 



7 l/3 



From the above VA equations for the 2D case, one can derive the value of A^ c 
for small values of e ^ as: N c — 1tt(1 — 8eexp(— 2)). For e — 0.2 we obtain 
N c — 9.845 which is in reasonable agreement with the value 10.8 obtained from 
numerical simulations of the 2D case. We remark that for small values of e the 
soliton is very extended in space and resembles a Bloch wave modulated with an 
envelope. In this case an effective mass approximation may be appropriate which 
allows to replace the GPE field equation with a nonlinear Schrodinger equation 
with effective mass m* and nonlinearity (3. In the 2D case we have 

iu t + a(u xx + u yy ) + (3\u\ 2 u — 0, (18) 

where a = m*/m,m<-^ = (d 2 E/dk 2 ) k=Q and (3 = {2tt/L 2 ) f d 2 r\<P h0 \ 4 (a 
similar equation can be written also for the 3D case) . In this approximation the 
norm of the gap Townes soliton can be evaluated as N — -^-g-. For deep optical 

lattices (e > 5) a can be approximated as aw e 1 / 4 and the norm N w Jijip ■ 



4 Numerical study of gap-Townes solitons and 
delocalizing transitions 

In this section we investigate the existence curve of gap-Townes solitons in the 
(N, e) plane. As mentioned before, this curve coincides with the delocalizing 
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transition curve which separates stable localized solutions from decaying ones. 
Its knowledge is therefore important for experimental investigations of multidi- 
mensional solitons. To this regard we remark that for the observation of multi- 
dimensional BEC solitons parameters should be chosen between the delocalizing 
and the collapsing curves. In the following we investigate the delocalizing curve 
by means of direct numerical integrations of the 2D and 3D GPE with periodic 
potential and by the corresponding ID NLS systems with quintic and septic non- 
linearities discussed above, respectively. The existence of gap-Townes solitons is 
then shown by direct numerical simulations of the multidimensional GPE using 
as initial conditions the above mentioned product states, which are found by 
solving the ID GPE with higher nonlinearities by means of the self-consistent 
method described in [20]. In Fig. [4] we depict the existence curve of gap-Townes 
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Fig. 4. Left panel: Existence curve of gap-Townes solitons of the 2D GPE for 
7 = 1. The continuous line represents the delocalizing curve as obtained from 
the ID quintic NLS approximation. The dots joined by thin dotted line represent 
numerical results obtained from direct integrations of the 2D GPE. Right panel: 
Existence curve of gap-Townes solitons of the 3D GPE for 7=1. The line repre- 
sents the delocalizing curve as obtained from the ID quintic NLS approximation, 
while the dots represent numerical results obtained from direct integrations of 
the 3D GPE. 



solitons of the 2D GPE for 7 = 1 as obtained from numerical integrations of the 
2D GPE. The corresponding curve obtained from the ID quintic NLS approx- 
imation by means of a self-consistent approach is also shown. We see that for 
e > 1 the ID quintic NLS curve agrees very well with that of the 2D GPE, the 
deviations becoming evident only for strengths of the OL which are less than 
£ w 1 (one recoil energy). This fact can be easily understood from the observa- 
tion that for a fixed value of e there is one value of N for which the gap Townes 
soliton exists and that by increasing e the corresponding value of N decreases. 
This implies that for gap-Townes solitons in a strong OL the nonlinear interac- 
tion is effectively small, due to the potential barriers which prevent tunneling 
of matter into adjacent wells. On the contrary, in a shallow optical lattice the 
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matter can easily tunnel through the barriers and the effective (attractive) non- 
linearity can be larger. Since the separability ansatz used to link the 2D GPE to 
the ID quintic NLS works well when the nonlinearity is small, it is clear that a 
discrepancy can arise at small values of e. The fact that the quintic NLS equation 
deviates from the 2D GPE only for e < 1, however, makes the mapping between 
these two equations very convenient for practical calculations. A similar situa- 




y y 

Fig. 5. Time evolution of the y-section of a 3D gap-Townes soliton obtained 
for e = 2 with slightly undercritical (left panel) and overcritical (right panel) 
number of atoms, corresponding to points a and b in FiglH respectively. 

tion seems to be true also for the 3D case. In the right panel of Fig. [4] we depict 
the existence curve of gap-Townes solitons of the 3D GPE for 7 = 1 as obtained 
from the ID septic NLS approximation by means of a self-consistent approach. 
Due to the long computational times required in 3D simulations we have pre- 
sented verification for only few points of the curve. A more complete analysis 
will require further investigation. In particular, the prediction of our analysis are 
checked by means of numerical integrations of the 3D GPE in correspondence of 
the two filled dots depicted in the right panel of Fig. 0] at e = 2, just above and 
below the delocalizing curve obtained from the septic ID NLS equation. In Fig. 
[5] we show the time evolution of the y-section of a 3D gap-Townes soliton ob- 
tained for e = 2 in correspondence of these slightly undercritical (left panel) and 
overcritical (right panel) points. We see that, in analogy to what we observed in 
2D case, the initial condition with an undercritical number of atoms leads to the 
complete delocalization of the state, while in the overcritical case a gap soliton 
state which remains stable over long time is formed. It is appropriate to mention 
that the exact profile of the 3D gap-Townes soliton, which separates these two 
behaviors is more difficult to obtain with the separability ansatz than for the 



12 Mario Salerno 1 , F. Kh. Abdullaev 2 , and B. B. Baizakov 2 



previous 2D case, since in the 3D case this ansatz appears to be less accurate. 
The signature of the gap-Townes soliton state, however, is very clear as one can 
see from the early stages of the time evolution depicted in Fig. [5] (notice from the 
left panel that the undercritical state remains stable up to a time t — 20 before 
starting to decay). This behavior strongly suggests the existence of gap-Townes 
solitons also in the 3D case. 

We now address to another important property of the gap-Townes soliton 
which shows up in the transition to a gap soliton state for overcritical number 
of atoms. As it was already mentioned, even a slightly overcritical norm gives 
rise to a rapid shrinking of the gap-Townes soliton. Although at initial stage of 
evolution this behavior is similar to the collapse of ordinary Townes soliton, the 
final state is different. Specifically, gap-Townes soliton approaches the (final) gap 
state via long-lasting oscillations, as shown in Fig. [6l Each "reflection" from the 
broad state is accompanied by emission of linear waves, which can be viewed 
as a tunneling of matter from the localized mode into the extended one. This 
process also contributes to the damping of the oscillations. We remark that these 
oscillations, having a very regular behavior, could be used to detect the existence 
of gap- Town solitons in a real experiment. 




2000 



Fig. 6. Oscillations of the gap-Townes soliton when it collapses into the final 
gap state, according to Eq. |T]) with initial state taken as the product solution 
constructed from two quintic NLS with a total norm increased by one percent. 
Shown is the y — cross section of the 2D profile. 



5 Conclusions 

In this paper we have shown the existence of gap-Townes solitons in the multidi- 
mensional (2D and 3D) Gross-Pitaevskii equation with a periodic potential and 
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discussed its link with the phenomenon of delocalizing transition. These solutions 
have the peculiarity of being unstable under small fluctuations of the number 
of atoms and separate localized (soliton like) states from extended (Bloch like) 
ones. The existence curve in the parameter space of this particular solution is 
very useful since it provides the lower threshold for the existence of localized 
states. The gap-Townes solitons discussed in this paper are a natural general- 
ization of the Townes solitons of the nonlinear Schrodinger equation (without 
periodic potential). The existence of Townes solitons in a nonlinear glass sam- 
ple modelled by the NLS equation was experimentally demonstrated in Ref. 
[21] . The fact that the transition from a gap-Townes soliton to a gap soliton 
is always accompanied by regular oscillations gives the possibility to indirectly 
observe the multidimensional gap-Townes solitons discussed in this paper in real 
experiments. In particular we expect these solitons to be observed both in mul- 
tidimensional BECs in OL and in nonlinear optics systems, including 2D and 
3D photonic crystals and arrays of nonlinear optical waveguides. 
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